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Abstract 

In this paper, the problem of terahertz pulsed imaging and reconstruction is 
addressed. It is assumed that an incomplete (subsampled) three dimensional 
THz data set has been acquired and the aim is to recover all missing samples. A 
sparsity-inducing approach is proposed for this purpose. First, a simple interpo¬ 
lation is applied to incomplete noisy data. Then, we propose a spatio-temporal 
dictionary learning method to obtain an appropriate sparse representation of 
data based on a joint sparse recovery algorithm. Then, using the sparse coeffi¬ 
cients and the learned dictionary, the 3D data is effectively denoised by minimiz¬ 
ing a simple cost function. We consider two types of terahertz data to evaluate 
the performance of the proposed approach; THz data acquired for a model 
sample with clear layered structures (e.g., a T-shape plastic sheet buried in a 
polythene pellet), and pharmaceutical tablet data (with low spatial resolution). 
The achieved signal-to-noise-ratio for reconstruction of T-shape data, from only 
5% observation was 19 dB. Moreover, the accuracies of obtained thickness and 
depth measurements for pharmaceutical tablet data after reconstruction from 
10% observation were 98.8%, and 99.9%, respectively. These results, along with 
chemical mapping analysis, presented at the end of this paper, confirm the ac¬ 
curacy of the proposed method. 
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1. Introduction 


Terahertz pulsed imaging (TPI) systems hold great potential in many appli¬ 
cations such as medical diagnosis of human tissue, the detection and chemical 
mapping of illicit drugs and explosives, and pharmaceutical tablet inspection 
[lj, M In a typical TPI measurement, the THz waveform for each pixel is 
recorded as a function of optical time delay. This provides a three-dimensional 
(3D) data set where two axes describe an object’s spatial content, and the third 
one represents the (time) depth information. In TPI, the transient electric field, 
rather than the radiation intensity, is measured. Thus, by applying the one di¬ 
mensional Fourier transform along the time axis we can get the magnitude and 
phase information of THz spectral data for each pixel. Despite these advantages, 
most existing TPI systems suffer from slow scanning speed and high implemen¬ 
tation cost due to pixel-by-pixel raster scan mechanism. One approach to speed 
up the acquisition process is to reduce the number of samples (measurements) 
and yet preserve the reconstruction quality. Recently, compressive sampling 
(CS) 0,0 has been widely used for this purpose. This theory essentially works 
based on random projections of input signals and sparse representation 0 dur¬ 
ing reconstruction. It holds great potential for sampling rates reduction, imag¬ 
ing time, power consumption and computational complexity. However, the main 
challenge in utilizing CS is implementation costs and specifically designing the 
sampling operator which requires costly and sophisticated hardware modules. 
Most recent researches on CS-THz has been focused to address this issue. For 
instance, in yl the authors propose a spining disk as a plausible sampling oper¬ 
ator for high-speed compressive acquisition. In 0, a block-based CS method is 
proposed, and in 0 a single pixel camera based on Bernoulli random matrix is 
presented. Recently, an improved version of 0, called adaptive CS, is reported 
10J. In this method, additional measurement points are adaptively added at the 
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regions prone to degradation with the aim of improving reconstruction quality. 

In this paper, we consider a very simple acquisition scenario, i.e. incomplete 
(subsampled) data, which can be easily realized without too much of hardware 
burden. Instead, we aim to utilize advanced reconstruction techniques to re¬ 
trieve the original data samples with least possible quality degradation. In the 
sequel, more details are explained. 

It should be noted that the spatial and temporal (or spectal) features in 3D 
TPI data have intrinsic geometrical structures. If these structures are modelled 
properly, it may lead to new TPI imaging systems with faster acquisition speed 
and more accurate data analysis. This aspect is an open issue and has not been 
specifically studied in the THz literature. Over the past few years, there have 
been increased interests in the study of sparse representation, in which a signal 
is characterized by a few non-zero coefficients in a certain transform domain. Up 
until now, sparse modelling has found applications in many image processing- 
related tasks such as acquistion, denoising, inpainting and super-resolution 111 . 


3Q,Q,Q. 


Most of existing works are limited to a deterministic sparsifying 
transform, such as the Fourier, the wavelet and the curvelet etc. More recently, 
it has been shown that a signal can be represented with fewer coefficients over 
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a learned dictionary from a number of training samples (see m and references 
therein). 

In this paper, we aim to investigate dictionary learning from incomplete and 
noisy 3D TPI data set. In particular, the missing data are first estimated from 
random subsets using Bicubic interpolation. Incomplete TPI data set can occur 
in two different scenarios, in which both require retrieving the missing samples; 
i) result of data transfer from a sender to receiver, and it) subsampled THz data: 
a simplified case of CS where merely some samples are missing (this obviously 
requires much simpler sampling operator than original CS). As a novel tool 
to address the sample recovery problem, 3D dictionary learning is used here. 
Spatial-temporal dictionary learning is applied to the 3D data set by exploiting 
the joint sparse model. Finally, denoising is performed based on the learned 
dictionary through convex optimization. Experimental results show that even 
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with 5% to 20% subsampled data, one can still get reliable spatial, structural 
and spectral information of the object. To the best of our knowledge, this is 
the first work demonstrating the advantages of 3D dictionary learning for THz 
data. These results can be exploited to speed up TPI measurement process 
substantially by reducing the quantity of data. 

This paper is organized as follows. In the next section we describe the 
proposed method including mathematical expression of the model, dictionary 
learning and joint sparse recovery. Section [3] is devoted to demonstrating its 
applicability to the experimental results. Finally, the conclusion is drawn in 
section [I] 


2. Proposed method 

2.1. Problem formulation 


The first step for recovering the original data from incomplete (subsampled) 
observations is inpainting. This can be simply achieved using a preliminary and 
fast approach such as Bicubic interpolation [8]. The main drawback is that the 
results of such techniques are of low quality especially in noisy scenarios. Other 
inpainting methods which rely on fixed transforms (e.g. 3D wavelet) [17] also do 
not lead to acceptable results as we will see later in our experimental results. A 
recent family of methods which try to build up an adaptive dictionary directly 
over the incomplete data has shown improved results for 2D natural images when 
the observation rate is above 25% [l^l ■ Efficient extension of these techniques for 
THz data, especially in severely subsampled noisy data (< 25%) is our objective. 
To take advantages of existing techniques for both incomplete and noisy data 
we propose a new model shown in Fig. [Q If we denote the original clean 
THz datacube by X £ ^N x xN y xB anc j incomplete noisy datacube by Y 
(of the same size), the estimated low-quality noisy datacube can be represented 
by Y (of the same size). Our aim is to denoise Y using a dictionary which is 
specifically designed to exploit both spatial and temporal correlations existing 
within THz data. 
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Figure 1: Block diagram of the proposed approach. 



Figure 2: Geometric illustration of block extraction procedure from the THz datacube. 


2.2. Dictionary learning 

Since the THz data is huge, one cannot learn a dictionary directly over the 
entire datacube. More importantly, 3D nature of THz data should be taken 
into account when constituting the training samples. Assume that y i is the 
i -th training vector which is obtained through partitioning Y into N small 3D 
blocks of size n x xn y x b , where, {n x , n y ) and b are the block sizes in spatial and 
temporal dimensions, respectively. These blocks may have spatial and temporal 
overlaps and are always converted to the column-vector of length r = n x n y b. 
Fig. [2] illustrates the geometrical view of the proposed block extraction. In 
addition to capturing spatio-temporal structure in the proposed scheme, we have 
the flexibility to define the contributions of spatial and temporal dimensions by 
adjusting ( n x ,n y ) and 6, respectively. 

The next step is to find the dictionary D £ R rxfe and sparse coefficients, 
denoted by {s J }^ =1 , from the training set {yi}^. In spite of most traditional 
dictionary learning methods in which all training vectors are treated indepen¬ 
dently, we devote a different approach. In order to effectively exploit the 3D 
spatio-temporal structure of THz data we cumulate every l training vectors into 
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a matrix f lj £l rxl and solve the following problem: 

j = D S j + E j , (1) 

where j denotes the j- th subset, D £ R rxfc is the dictionary, S j £ M. kxl is the 
corresponding sparse coefficients matrix, and Ej £W rxl represents the decom¬ 
position error. This is a joint sparse model which simultaneously deals with the 
training vectors in f lj. There are different ways to choose the training subset 
ftj. For instance, one can cluster the 3D blocks and then group them based on 
the amount of similarity within the voxels of these blocks. For simplicity, and yet 
preserving the sparsity structure of neighboring blocks, we use raster-scanning 
with full spatial and temporal overlap from left-to-right and top-to-bottom as 
shown in Fig. [2] We have observed encouraging results using this scheme, but 
the performance of more sophisticated grouping techniques can be further in¬ 
vestigated in the future. Here, we group all the blocks into \N/t] subsets of size 
l: 

{yi, • ■ • yz}, {y/+i, ■ ■ • y2 1 }, ■ ■ ■ { yn - i + i , ■ ■ ■ y^}, 

Si v ✓ S y ✓ S v ✓ 

S7i n riV / n 

where [•] is the ceiling operator. Correspondingly, S ? ! s are defined as k x l 
matrices as follows: 


{ s l, • • • S(}, {s; + i, . . . S 2 /}, . . . {sjV-J+1, • ■ • S N }, 


wn 


and the aim is to solve the following minimization problem subject to sparsity 


of Sj\ 


Vj min \\flj — DS 


o\\ 2 f 


( 2 ) 


Most typical dictionary learning approaches work based on “alternating min¬ 
imization” [hi a [ 20 ! . These methods mainly contain two major steps, co¬ 
efficient update and dictionary update, which are alternately executed until 
reaching local minima for P]|. 

In order to jointly update the sparse coefficients we choose MMV (multiple 
measurement vector) sparse coding framework. In this framework, as opposed 
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Figure 3: MMV model 


to SMV (single measurement vector), several sparse vectors are simultaneously 
recovered. The main assumption in MMV is that the sparse vectors admit a 
common sparsity pattern, i.e., the locations of non-zeros are the same for all 
vectors. Fig. [3] is a diagram illustrating the MMV model. There are several 
algorithms in the literature to solve 0 with respect to S , . Most of these 
methods are extended from SMV-based approaches [ 21 I, [ 22 , 23 L We have found 
that SOMlj* (Simultaneous Orthogonal Matching Pursuit) 22, 231 works well 


for reconstruction of THz data and hence we utilized it for coefficients update. 

Next, for the dictionary update (solving 0 with respect to D), we use a 
well-established method called K-SVD Q. K-SVlJf| is a generalization of K- 
means clustering and updates the dictionary in a column-by-column scheme by 
using singular value decomposition (SVD) (full details in }16l|). 


2.3. THz data reconstruction 

After finding D and {s 1; we need to estimate the THz data. To do this, 
we define the following minimization problem: 

N N 

min A ||y - ^ ll R * x “ Ds i II2 + P H R * X ~ m ‘ll2 • ( 3 ) 

i =1 i=l 


1 Available in SPAM: http://spams-devel.gforge. inria.fr/doc/html/doc_spams002.html 
2 KSVD-Box: http: //www. cs. technion. ac . il/~ronrubin/software .html 
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where y = vec(Y) and x = vec(X), both with length p = N x N y B , are vectorized 
representations of datacubes Y (noisy) and X (clean), respectively!^ Also, y 7 ; = 
R;y and x, = R^x, where R^ is a huge binary matrix of size r x p. R, has only 
one non-zero (i.e. 1) in each row. As a block extraction operator, R, extracts 
the voxels belonging to i-th block. 

In J3]), the leftmost term is the error between noisy and clean data, the 
middle term is related to the sparse decompostion error, and the rightmost 
term is the smoothness constraint with regularization parameter /3. Moreover, 
m, = [ rrii , mj,...] 3 is the r x 1 mean vector of the i-th block (i.e. average of all 
elements in R^x). Problem (U) is convex in x and can be solved by zeroing its 
gradient, with respect to x, which finally leads to: 

x= ^AI+(l+/3)^RfR^ ( A y + XX (Dsi+ZIm^ . (4) 

Although the above expression seems computationally expensive at the first 
glance, it does not need to be directly calculated in practice. Instead, since the 
inverting term in (U) is diagonal, we obtain the estimated datacube in a voxel- 
wise fashion (similar to the strategy used in Q) . We can show the voxel-wise 
calculation of one voxel x (taking the corresponding voxels in y and m into 
account) using the following sets of equations: 


min \C(x) = X(x — y ) 2 + (x — ds ) 2 + /3(x — m) 2 ) 
dC 

—> — = 2\{x — y) + 2{x — ds) + 2/3(x — m) = 0 


1 


—> x = 


(A y + ds + /3m) 


(5) 


1 + /3 + A 

where C is the cost function to be minimized, and ds is the corresponding voxel 
in Ds. It is seen that the last equation above well matches with (U) which is 
the non-practical form of reconstructing datacube. 

The pseudo-code of the proposed method is given in Algorithm [Q 


3 Note that x and y should not be confused with 3D blocks where always have subscript 
indices thorough the paper. 
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Algorithm 1: Pseudo-code of the proposed method. 


Input: Incomplete noisy datacube Y . 
Output: Recovered datacube X. 

l begin 


2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 


Interpolation: 

Apply Bicubic interpolation on Y and yield Y; 

Dictionary learning: 

Extract the spatio-temporal blocks from Y and obtain {yi} i=1 ; 

Group the blocks and constitute 

Initialize the dictionary D with the DCT (discrete cosine transform) basis; 

repeat 

Apply SOMP to flj and yield S j for all j = 1. ... 

Update D using K-SVD dictionary update |lg |: 
until stopping criterion is met; 

THz data reconstruction : 

Calculate Q to yield x and then reshape it to X; 


14 end 


3. Experimental results 

We present the results of applying the proposed method to two different 
types of THz data. The first dataset, which we call it “T-shape”, has size 
200 x 200 x 512, and acquired across an area of 20mmx22mm using a TPIscan- 
1000 system (TeraView Ltd, Cambridge, U.K.), covering a spectral range from 
0.1 to 3.5 THz. The sample used is a polythene pellet of a diameter of 25 mm. 
Inside the pellet there is a T-shaped plastic sheet which locates approximately 
0.2 mm below the sample surface. For this dataset, the structural information 
such as thickness and depth are of interest. The second data is a low spatial 
resolution data called “Tablet data” acquired from pharmaceutical tablets. Two 
of such sets (called LA and TP) have size 42 x 42 x 512, and other seven sets 
have size 49 x 49 x 609. Extraction of spectral information such as chemical 
mapping is crucial for this dataset. 
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(a) (b) 


SNR=16.23 dB SNR = 19.02 dB 




(c) (d) (e) 

Figure 4: T-shape data at one temporal band: (a) original image, (b) 10% incomplete noisy 
observations. Reconstruction results using (c) soft-thresholding in 3D wavelet domain, (d) 
Bicubic interpolation, and (e) proposed method. 

3.1. T-Shape data 

Following the proposed model, we first applied Bicubic interpolation to the 
noisy incomplete data, yielding Y. We considered input additive Gaussian noise, 
leading to input SNR in the range 15 ~ 20 dB, which is a reasonable range in 
practice. We then applied the proposed dictionary learning method to find the 
dictionary and corresponding sparse coefficients to reconstruct the entire dat- 
acube. Due to importance of structural information in T-shape data, we gave 
more contribution to the spatial domain than temporal domain for the dictio¬ 
nary learning. Hence, the experiments were conducted by setting the spatial 
size n x = n v = 8, and two different temporal sizes: 6=1 (spatial dictionary), 
and 6 = 4 (spatio-temporal dictionary). In any case we used full spatial and 
temporal overlaps. According to these settings, the obtained dictionaries were 
of size 64 x 256 and 256 x 256. A sample reconstructed image with 6 = 4, l = 10, 
A = 0.5, /3 = 0.1, and input SNR of 17 dB is given in Fig. [4] It is seen from 



SNR = 16.21 dB 



10 






Figure 5: T-shape data: SNR of reconstructed data at 10% observation rate versus input 
SNR. 


this figure that the proposed method has recovered and denoised the original 
image from only 10% of noisy samples with highest SNR among other methods. 
The proposed method was able to enhance the SNR of Bicubic interpolation 
results by up to 2.5 dB. Other methods, i.e. soft-thresholding after applying 
3D symlet4 wavelet transform [l7|, shows weaker performances compared to the 
proposed method. We also show in Fig. [5] the reconstructed SNR at the fixed 
observation rate of 10% versus different input SNRs. As seen from this figure, 
the proposed approach outperforms other methods. 

Next, performance of the proposed method against different values of b and 
l is investigated. For this purpose, SNRs of the recovered datacube for different 
observation rates are given in Table [T] From this table, the advantage of spatio- 
temporal dictionary can be observed by comparing the SNRs at b = 4, with 
those at b = 1. Also, the given SNRs in this table indicates that the proposed 
method is more effective than 3D wavelet transform. We empirically observed 
that l « 10 gives the best performance for T-shape dataset. 

Tabic [2] gives an insight about the computational complexity of the proposed 
method. In this table, the elapsed time of the dictionary learning stage for 
different selections of b and l is shown. It is clearly seen that large values of 
l can significantly reduce the computation time. This behavior supports the 
advantages of using joint sparse recovery which has already shown to improve 
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Table 1: T-shape data: recovery accuracy (SNR) in dB for different methods. The results are 
shown for several observation rates, and different selections of b and l. 


observation rate 

SNR (dB) 

b, l 

3D wavelet 

Bicubic interp. 

proposed method 

5% 

10.53 

16.89 

17.45 

1, 1 

17.69 

1, 10 

18.78 

4, 1 

18.96 

4, 10 

10% 

14.71 

17.68 

19.41 

1, 1 

19.46 

1, 10 

20.17 

4, 1 

20.94 

4, 10 

15% 

16.72 

17.87 

20.05 

1, 1 

20.33 

1, 10 

20.86 

4, 1 

22.25 

4, 10 

20% 

17.69 

17.96 

20.88 

1, 1 

21.03 

1, 10 

22.53 

4, 1 

23.15 

4, 10 


Table 2: T-shape data: The computation time of dictionary learning step in the proposed 
method for different selections of b and /. These times are in second and calculated per frame. 


o bservat ion 
rate 

b, l 

10% 

20% 

1,1 

110.3 

122.5 

1, 10 

5.0 

7.6 

4, 1 

93.3 

97.6 

4, 10 

4.8 

5.9 


the quality as well (in Table [T]). It can be further seen from Table [2] that 
increasing b also decrease the computation time, though not as dramatically as 
that for l. 

As a useful measure, we performed the thickness and depth evaluations be¬ 
fore and after reconstruction of T-shape object. These results give structural 
information about the object of interest. Fig. [6] (a) geometrically illustrates 
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(a) (b) 

Figure 6: T-shape data: (a) geometrical representation of T-shape object, and (b) sample 
time-domain waveform at one spatial location. 


Table 3: T-shape data: measured “meamtstandard deviation” of thickness and depth for 
different observation rates._ 



observation rate 

reconstructed 

original 

Thickness (mm) 

5% 

0.0766 ± 0.0063 

0.0750 ± 0.0029 

10% 

0.0741 ± 0.0063 

15% 

0.0756 ± 0.0056 

20% 

0.0748 ± 0.0054 

Depth (mm) 

5% 

0.1889 ±0.0079 

0.1902 ±0.0071 

10% 

0.1909 ±0.0087 

15% 

0.1911 ±0.0079 

20% 

0.1899 ±0.0082 


these two parameters. Thickness (t) and depth (d) calculations are practically 


Q. 


useful to identify inhomogeneities or defects in the object of interest [If. In order 
to find these parameters we refer to a sample temporal waveform at one spatial 
location shown in Fig. [3] (b). The obtained values for these two parameters 
are subject to appropriate scalings to be converted to millimeter (more details 
can be found in [l|). The calculated means and standard deviations of these 
parameters for different observation rates are given in Table [3] The results in 
this table are shown for both original and reconstructed data using the proposed 
method. It is found from Table [3] that the accuracy of reconstruction using the 
proposed method is very high, as the depth and thickness are approximately 
equivalent to the original data. 
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LA TP 



Figure 7: Tablet data: the reconstructed SNR versus observation rate for both LA and TP 
datasets at 20 dB input SNR. 


3.2. Tablet data 

The dictionary learning settings for Tablet data is different from those of 
T-Shape data. Due to very low spatial resolution of Tablet data, it is natu¬ 
ral to give a much higher contribution to the temporal dimension during the 
dictionary learning. This is inline with the fact that the temporal/spectral in¬ 
formation is more valuable than spatial information in pharmaceutical studies 
[1] 251. Therefore, we selected n x = n y = 2 and b = 128 with full spatial and 
temporal overlaps within the blocks. 

As a quantitative measure, we show in Fig. [T] the variations of reconstructed 
SNR of both LA and TP datasets versus different observation rates, where 
l = 10, k = 64, A = 0.5 and /3 = 0.2 were used for the proposed method. 
These values are chosen empirically and are manually tuned to achieve the best 
performance among other selections. We also observed that the achieved results 
are not sensitive to variations of A and [3 around the selected values. The 
dictionary columns size, i.e. fc, is set to 64 due to lower spatial resolution in 
Tablet data compared to T-shape data. The input SNR for this experiment was 
set to 20 dB. Similarly, we show the reconstructed SNR versus input SNR at 
20% observation rate. It is clearly seen that the proposed method is able to 
significantly improve the reconstructed SNR in both figures. 

In order to observe the robustness of the proposed approach, more sets of 
tablet data were included in our experiments and the average performance was 


14 

















LA TP 



5 10 15 20 25 30 5 10 15 20 25 30 

Input SNR (dB) Input SNR (dB) 


Figure 8: Tablet data: the reconstructed SNR versus input SNR for both LA and TP datasets 
at 20% observation rate. 


measured. Seven incomplete Tablet datasets at different observation rates were 
reconstructed by using well-known methods. Three input SNR levels were used 
in this experiment. The results are given in Table [4] It is seen that the achieved 
results using the proposed method is better than those obtained using conven¬ 
tional techniques. Moreover, the average output SNRs are compatible with what 
obtained in the previous experiment (Figures [7] and [5]). 

As mentioned before, chemical mapping is a useful tool to analyze and ob¬ 
serve the uniformity of the tablet Q. Chemical mapping can be obtained using 
the terahertz spectral information. This evaluation should be performed in 


Table 4: Tablet data: output SNRs (dB) at different observation rates and input SNRs (dB), 


averaged over seven datasets. 


Input SNR (dB) 

observation rate (%) 

Output SNR (dB) 

3D wavelet 

Bicubic interp. 

proposed method 


10 

7.91 

9.49 

14.71 

10 

20 

9.31 

11.14 

16.26 


30 

12.16 

13.75 

18.12 


10 

8.34 

12.40 

16.13 

20 

20 

10.67 

14.83 

17.22 


30 

12.24 

15.91 

19.49 


10 

9.81 

14.92 

18.23 

30 

20 

12.64 

16.78 

20.12 


30 

14.31 

18.63 

22.70 
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the frequency-domain, thus, we first take Fourier transform of the time-domain 
waveforms for both original and reconstructed data. Then, a reference (spec¬ 
tral) waveform denoted by u is selected from original LA dataset at one pixel 
location (normally from the middle of the tablet). After that, a spectral match¬ 
ing should be performed to generate the chemical map. We used the cosine 


correlation mapping (CCM) |26|, as used in 1], for this purpose. If we define a 
spectral frequency-domain vector by v, then the CCM can be calculated via: 

p 


cos 


<«) = £ 


^ p * Vp 


,.UL • 

p— 1 11 11 2 


( 6 ) 


where P is the number of spectral bands. Parameter 9 is the angle between 
v and u, and therefore, cos (9) is between 0 and 1. This evaluation should be 
performed for all pixels of the corresponding dataset. Smaller angle (larger 
cos (9)) means better match between reference spectra and the waveform under 
evaluation. The results of chemical mapping for our tablet data is shown in 
Fig. [9] We calculated CCM for four cases; TP with LA reference, LA with LA 
reference, TP with TP reference, and TP with LA reference. It is seen from 
Fig. [9] that both the reconstructed data have a close chemical map similar to 
the original ones, and that they have a smooth and uniform shape. 


4. Discussion and Conclusion 

In this paper, we addressed the THz data recovery from an incomplete noisy 
set of observations. The main advantage of starting with incomplete (subsam¬ 
pled) THz data is fast and inexpensive acquisition process. Inspired by recent 
developments in dictionary learning and sparse recovery frameworks, we pro¬ 
posed a spatio-temporal dictionary learning technique to find an efficient sparse 
domain for the THz data. The proposed method aimed at exploiting the exist¬ 
ing temporal and spatial correlations. It is worthy here to know the difference 
between the time-domain THz system and visible-light camera by noting their 
image formation details. In fact, even for a static object, one can still get 3D 
THz data, while for a visible-light camera, the 3D data leads to a video sequence 
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TP (left) and LA (right) with TP reference LA (left) and TP (right) with LA reference 



10 20 30 40 50 60 70 80 10 20 30 40 50 60 70 80 


(b) 

Figure 9: Tablet data: the chemical map for (a) original and (b) reconstructed data from 20% 
observation rate. 

with moving scene. To be more accurate, in terms of spatial correlation, THz 
signal exhibits some similar correlations as that obtained by a visible camera. 
However, in terms of temporal-correlation, it is significantly different from that 
of video sequences at visible spectrums. Note that the temporal THz data will 
either give about structural information like thickness and depth or provide 
some spectral information, e.g., the chemical component of an object, as pre¬ 
sented in the previous section. Extensive sets of experiments were conducted to 
support the effectiveness and accuracy of the proposed method. The thickness 
and depth calculations, and chemical mapping analysis for two different types 
of THz data, along with the achieved SNRs, confirmed the advantages of the 
proposed approach. In future, we are going to extend the proposed method 
for the complex scenario where a joint spatio-temporal and frequency domain 
dictionary learning can be achieved. 
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